home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / src / qr.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1997-07-10  |  3.4 KB  |  162 lines
  1. /*
  2.  
  3. Copyright (C) 1996 John W. Eaton
  4.  
  5. This file is part of Octave.
  6.  
  7. Octave is free software; you can redistribute it and/or modify it
  8. under the terms of the GNU General Public License as published by the
  9. Free Software Foundation; either version 2, or (at your option) any
  10. later version.
  11.  
  12. Octave is distributed in the hope that it will be useful, but WITHOUT
  13. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15. for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with Octave; see the file COPYING.  If not, write to the Free
  19. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  20.  
  21. */
  22.  
  23. #ifdef HAVE_CONFIG_H
  24. #include <config.h>
  25. #endif
  26.  
  27. #include "CmplxQR.h"
  28. #include "CmplxQRP.h"
  29. #include "dbleQR.h"
  30. #include "dbleQRP.h"
  31.  
  32. #include "defun-dld.h"
  33. #include "error.h"
  34. #include "gripes.h"
  35. #include "help.h"
  36. #include "oct-obj.h"
  37. #include "utils.h"
  38.  
  39. DEFUN_DLD (qr, args, nargout,
  40.   "[Q, R] = qr (X):      form Q unitary and R upper triangular such\n\
  41.                        that Q * R = X\n\
  42. \n\
  43. [Q, R] = qr (X, 0):    form the economy decomposition such that if X is\n\
  44.                        m by n then only the first n columns of Q are\n\
  45.                        computed.\n\
  46. \n\
  47. [Q, R, P] = qr (X):    form QRP factorization of X where\n\
  48.                        P is a permutation matrix such that\n\
  49.                        A * P = Q * R\n\
  50. \n\
  51. [Q, R, P] = qr (X, 0): form the economy decomposition with \n\
  52.                        permutation vector P such that Q * R = X (:, P)\n\
  53. \n\
  54. qr (X) alone returns the output of the LAPACK routine dgeqrf, such\n\
  55. that R = triu (qr (X))")
  56. {
  57.   octave_value_list retval;
  58.  
  59.   int nargin = args.length ();
  60.  
  61.   if (nargin != 1 && nargin != 2 || nargout > 3)
  62.     {
  63.       print_usage ("qr");
  64.       return retval;
  65.     }
  66.  
  67.   octave_value arg = args(0);
  68.  
  69.   int arg_is_empty = empty_arg ("qr", arg.rows (), arg.columns ());
  70.  
  71.   if (arg_is_empty < 0)
  72.     return retval;
  73.   else if (arg_is_empty > 0)
  74.     return octave_value_list (3, Matrix ());
  75.  
  76.   QR::type type = (nargout == 0 || nargout == 1) ? QR::raw
  77.     : (nargin == 2 ? QR::economy : QR::std);
  78.  
  79.   if (arg.is_real_type ())
  80.     {
  81.       Matrix m = arg.matrix_value ();
  82.  
  83.       if (! error_state)
  84.     {
  85.       switch (nargout)
  86.         {
  87.         case 0:
  88.         case 1:
  89.           {
  90.         QR fact (m, type);
  91.         retval(0) = fact.R ();
  92.           }
  93.           break;
  94.  
  95.         case 2:
  96.           {
  97.         QR fact (m, type);
  98.         retval(1) = fact.R ();
  99.         retval(0) = fact.Q ();
  100.           }
  101.           break;
  102.  
  103.         default:
  104.           {
  105.         QRP fact (m, type);
  106.         retval(2) = fact.P ();
  107.         retval(1) = fact.R ();
  108.         retval(0) = fact.Q ();
  109.           }
  110.           break;
  111.         }
  112.     }
  113.     }
  114.   else if (arg.is_complex_type ())
  115.     {
  116.       ComplexMatrix m = arg.complex_matrix_value ();
  117.  
  118.       if (! error_state)
  119.     {
  120.       switch (nargout)
  121.         {
  122.         case 0:
  123.         case 1:
  124.           {
  125.         ComplexQR fact (m, type);
  126.         retval(0) = fact.R ();
  127.           }
  128.           break;
  129.  
  130.         case 2:
  131.           {
  132.         ComplexQR fact (m, type);
  133.         retval(1) = fact.R ();
  134.         retval(0) = fact.Q ();
  135.           }
  136.           break;
  137.  
  138.         default:
  139.           {
  140.         ComplexQRP fact (m, type);
  141.         retval(2) = fact.P ();
  142.         retval(1) = fact.R ();
  143.         retval(0) = fact.Q ();
  144.           }
  145.           break;
  146.         }
  147.     }
  148.     }
  149.   else
  150.     {
  151.       gripe_wrong_type_arg ("qr", arg);
  152.     }
  153.  
  154.   return retval;
  155. }
  156.  
  157. /*
  158. ;;; Local Variables: ***
  159. ;;; mode: C++ ***
  160. ;;; End: ***
  161. */
  162.